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The Diffusion Monte Carlo (DMC) method is applied to compute the ground state energies of the water 
monomer and dimer and their D 2 O isotopomers using MB-pol; the most recent and most accurate ab inito- 
based potential energy surface (PES). MB-pol has already demonstrated excellent agreement with high level 
electronic structure data, as well as agreement with some experimental, spectroscopic, and thermodynamic 
data. Here, the DMC binding energies of (H20)2 and (D20)2 agree with the corresponding values obtained 
from velocity map imaging within, respectively, 0.01 and 0.02 kcal/mol. This work adds two more valuable 
data points that highlight the accuracy of the MB-pol PES. 


The development of a full-dimensional potential energy 
surface (PES) for a many-body system that extends to 
progressively larger cluster sizes and, at the same time, 
is computationally feasible has been a longstanding issue 
in electronic structure theorjl^. Many different empir¬ 
ical water PESs have been parametrized ranging from 
fully coarse-grained to flexible atomistic models that in¬ 
clude polarization effects and charge transfer (see, e.g., 
Refs.[2Hni). Nevertheless, empirical PESs are generally in¬ 
adequate for capturing the behavior of water for a wide 
array of cluster sizes (i.e., from small clusters to the bulk 
liquid) and thus, have generated ample motivation for the 
construction of ab initio and ab inito-hased surfaces with 
the latter having a foundation in the many-body expan¬ 
sion of the interaction^i^. Several notable PESs belong¬ 
ing to this family are ppp jHl CC-pofp^, the HBBO-2 
series for the water dimeJi^l^i^, WHBl^iS, and HBB2- 
po }i7 | i8 | Along the same vein, the MB-pol PE^^^HUI has 
most recently emerged as an ab inito-hased water surface 
rigorously derived from the many-body expansion of the 
interaction energy and expressed in terms of explicit one-, 
two-, and three-body contributions, with all higher-order 
terms being represented by (classical) many-body induc¬ 
tion within a modified version of the polarization model 
originally employed by the TTM4-F potentiaP. Simi¬ 
larly to WHBB and HBB2-pol, the MB-pol one-, two-, 
and three-body terms were obtained from fits to large 
sets of CCSD(T) monomer, dimer, and trimer energies 
calculated in the complete basis set limit. 

Paesani and co-workers (see Refs. [MSI]) have demon¬ 
strated that MB-pol effectively reproduces high-level ab 
initio data for: (a) stationary point energies on the 
(H20)2 and (H20)3 PESs; (b) PESs for (H20)2 and 
(H 20)3 interaction energies plotted versus 0-0 distance 
and 0-0-0 angle, respectively; and (c) (H 20)4 inter¬ 
action energies and (H20)4_6 isomer energies. Like¬ 
wise, results from MB-pol were shown to exhibit good 
agreement with experimental data for (H 20)2 vibration- 
rotation tunneling splittingji^, as well as the structural, 
thermodynamic and dynamical properties of bulk water 
at a fully quantum-mechanical level^. Both infrared and 
Raman spectra calculated from centroid molecular dy¬ 
namics simulations with the MB-pol potential were found 


to be in good agree ment with the corresponding experi¬ 
mental results. 1221221 

In this paper, we use the diffusion Monte Carlo (DMC) 
method to compute the MB-pol (H20)2 and (D20)2 
binding energies and co mpar e them to the correspond¬ 
ing experimental values^SESl. Note here that Ref. Uni 
reported a DMC result for the binding energy of the 
HBB2 water dimer, (H20)2, that was later confirmed 

experimentallji2il. 

DMC uses a population of Nw random walkers that 
sample the configuration space bounded by a PES and 
collectively represent the wavefunction of the many-body 
system at projection time r for each time step of length 
A 7 ^ 2 lIIZl_ At sufficiently long r the distribution of ran¬ 
dom walkers becomes stationary and the instantaneous 
energy E^ei{T) fluctuates about its average value (Eref)- 
In the limit of r — >• 00 , At — >• 0, and —>• 00 , this 

distribution converges to the ground state wavefunction 
with (Ej-ef) = Eq, the ground state energy. In a re¬ 
cent papep2SI, we have undertaken a thorough analysis 
of the behavior and extent of the bias (systematic error) 
arising from the finite time step At, as well as the bias 
caused by a finite random walker population N\y for the 
q-TIP4P/El^ water monomer, dimer and hexamer. The 
time step bias in the estimate of Eq vanishes slowly, and 
for the water monomer at At = 10.0 au is 0.015 kcal/mol. 
However, this bias cancels nearly completely when the 
ground state energy differences, such as the binding en¬ 
ergy, 

Do ■= 2Eii^o - £’(H20)2! (1) 

or the isotope shift, 

(MJo :=Do(D20)-Do(H20), (2) 

are computed. (Note that the same value of At = 10.0 
au was used in the DMC calculations for the WHBB 
water hexameil^.) The bias in Nw also cancels for the 
isotope shift (SDq), but for the binding energy {Dq) it 
does not cancel. This is because the DMC ground state 
energy estimate for the monomer converges much faster 
with respect to Nyy than that for the dimer. Thus, in 
order to obtain an accurate estimate of the asymptotic 
value of Do at Nw —>■ 00 , one needs to perform a series of 










2 


calculations using several different values of Ny/- It was 
concluded that 0.002 kcal/mol accuracy for the dimer 
binding energy Dq could be achieved using ~ 2.0 x 
10^. 

In this work, we compute the energies for the MB-pol 
water monomer and dimer adhering to the DMC protocol 
of Ref. [55] and performing studies for the time step with 
Ar = 5, 10, and 15 au and for the random walker popula¬ 
tion with Nw = 1.6 X10^, 3.6 x 10^, and 1.96 x 10^. All of 
the DMC simulations were run to a maximum projection 
time of 2.0 x 10® au. In the time step studies, the walker 
population was fixed at Ny/ = 1.96 x 10'^, while a time 
step of At = 10.0 au was used in every walker number 
study. The asymptotic behavior of the DMC estimate of 
the ground state energy, Eq{At), in the Ar —>■ 0 limit 
can be approximated well by a polynomial with a vanish¬ 
ing derivative at Ar = 0. The latter is because the error 
caused by the split-operator approximation implemented 
in the DMC method is quadratic in At. Consequently, 
since three different values of the time step At have been 
used, we consider the cubic interpolation: 

Eo{At)^A + BAt'^ + CAt^, (3) 

using three fitting parameters, with A giving the ground 
state energy estimate. 

Here, we also comment that our walker number stud¬ 
ies employ larger values of Ny/ than might seem to be 
needed, as including more walkers in the population si¬ 
multaneously reduces the statistical uncertainty as well 
as the systematic erroil^. 

Fig. [T] and Fig. show the bias in At and Ny/ for 
the ground state energy estimates from DMC (Ag) as 
functions of the time step and the inverse walker popula¬ 
tion 1/Nw, respectively, as well as the associated bind¬ 
ing energies (Dq). The curves follow a pattern similar 
to that reported in Ref. [28] for the q-TIP4P/F dimer 
and monomer. A clear bias does exist for the ground 
state energies (Eq) in both At and Nyr, but for the re¬ 
ported range of time steps and walker numbers, the en¬ 
ergy change is in the second or third positions beyond 
the decimal. Moreover, the binding energy bias in At 
is virtually negligible due to error cancellations. This is 
because the largest contribution to the time step error 
comes from the intramolecular degrees of freedom, which 
are essentially the same regardless of whether the sys¬ 
tem is comprised of a single water molecule or multiple, 
interacting water molecules. Because the intramolecular 
modes are nearly invariant to changes in the number and 
spatial orientation of the water molecules in clusters, the 
behavior and extent of the Eg bias curves in At for the 
H 2 O and D 2 O monomer and dimer are very similar (see 
the top four panels of Fig. [^, which leads to substan¬ 
tial elimination of the bias such that the binding energy 
displays only a weak dependence on At. On the other 
hand, the bias in 1/Nw persists for the binding energy 
even after the energy differences are taken. In this case, 
the Eg monomer bias in Ny- disappears much faster than 
that of the dimer (note the difference in scales between 


the monomer (top) and dimer (middle) panels of Fig. 
[^, thereby causing imperfect cancellation of error and 
the appearance of a strong residual bias in the binding 
energy. 

Table. |l] shows the DMC binding energies extrapo¬ 
lated to the Ny —>■ 00 limit for three different PESs ex¬ 
amined in this work: q-TIP4P/I^, TTM3/I®, and MB- 
poiEMU] The statistical errors for the binding energies 
with the largest walker population used in this study of 
Ny = 1.96 X 10"^ are all known to be on the order of 
10“® kcal/mol, but these numbers are not displayed in 
the table because extrapolation of the error bar values 
has proven to be unreliable. Here, the same DMC proce¬ 
dure was used to calculate the ground state energies for 
the TTM3/F H 2 O and D 2 O monomer and dimer as those 
with q-TIP4P/F and MB-pol. Additionally, we also pro¬ 
vide the (H20)2 TTM3/F and HBB2 binding energies as 
reported in Ref. [151 The Dg values for the q-TIP4P/F 
dimers are ^ 1.4 kcal/mol larger than the experimen¬ 
tal bin ding e nergies obtained from velocity map imaging 
(Expt.JpSES, This discrepancy is not surprising consid¬ 
ering that the q-TIP4P/F PES was not designed to rep¬ 
resent the microscopic properties of small water clusters 
accurately. Our TTM3-F Dg value for (H20)2 agrees well 
with that reported by Bowman and co-workergl^ with a 
difference of only 0.06 kcal/mol, but our results indicate 
that the binding energy for this PES is still off from the 
experimental values by 0.62 and 0.53 kcal/mol for (H20)2 
and (D 20 ) 2 , respectively. 

Notably, the (H20)2 DMC binding energy, also from 
the ab initio-hased two-body PES, HBBilS, (with values 
of At and Ny similar to those used in the present work) 
is highly accurate upon comparison to the experimental 
results. At the same time, we highlight that the MB-pol 
PES likewise shows nearly perfect agreement with the 
velocity map imaging Dg values for both MB-pol dimers, 
(H20)2 and (D20)2. For the MB-pol PES, the binding 
energies differ from the experimental results by only 0.01 
kcal/mol for (H20)2 (as is also the case for the HBB2 
dimer PES) and 0.02 kcal/mol for (D20)2. Such an excel¬ 
lent correspondence adds two more valuable data points 
which underscore a favorable assessment for the accuracy 
of the MB-pol PES. 

Therefore, efforts are currently underway to establish 
the true ground state energy and wavefunction of the 
MB-pol water hexamer using DMC. Additionally, in ac¬ 
cordance with Ref. [551 the ground state energies for the 
isomers of the MB-pol hexamer will be computed. 
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FIG. 1. DMC energies for the MB-pol water monomer and dimer as a function of the time step At. The time steps used in 
this study were At — 5.0, 10.0, and 15.0 au. The walker population was fixed at Nw = 1-96 x 10"^, and the total projection 
time for all runs was 2.0 x 10® au. Data points were interpolated using a cubic polynomial fit with a vanishing derivative at 
At — > 0. Top left: H 2 O monomer energies. Top right: D 2 O monomer energies. Middle left: H 2 O dimer energies. Middle right: 
D 2 O dimer energies. Bottom left: H 2 O dimer binding energy. Bottom right: D 2 O dimer binding energy. 






TABLE I. DMC binding energy Do for the H 2 O and D 2 O dimer with different PESs. The DMC binding energies computed 
in this work for a time step of At = 10.0 au were extrapolated to the Nw —> 00 limit. Bindi ng en ergies marked with the “a” 
superscript were reported in ref.^. The experimental binding energies (Expt.) are from refs.^^^^. The error bar magnitudes 
for this work are on the order of 10“®. All energies are reported in kcal/mol. 


Do 


Structure 

q-TIP4P/F 

TTM3/F 

MB-pol 

HBB2 

Expt. 

(H20)2 

4.53 

3.78, 3.84 ±0.07“ 

3.15 

3.15 ±0.01“ 

3.16 ±0.03 

(D20)2 

4.97 

4.09 

3.54 

3.56 ±0.03 
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FIG. 2. Same as Fig. [^but as a function of reciprocal walker number 1/Nw- The walker populations used in this study were 
Nw = 1.6 X 10^, 3.6 X 10^, and 1.96 x 10'*. The time step was fixed at Ar = 10.0 au, and no special interpolation technique 
was employed. 
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